function EDD = emaxdrawdown(Mu, Sigma, T)
%EMAXDRAWDOWN Compute expected maximum drawdown for a Brownian motion.
%
%	EDD = emaxdrawdown(Mu, Sigma, T);
%
% Inputs:
%	Mu - The drift term of a Brownian motion with drift.
%	Sigma - The diffusion term of a Brownian motion with drift.
%	T - A time period of interest or a vector of times.
%
% Outputs:
%	EDD - Expected maximum drawdown for each time period in T.
%
% Notes:
%	If a geometric Brownian motion with stochastic differential equation
%		dS(t) = Mu0 * S(t) * dt + Sigma0 * S(t) * dW(t) ,
%	convert to the form here by Ito's Lemma with X(t) = log(S(t)) such that
%		Mu = M0 - 0.5 * Sigma0^2
%		Sigma = Sigma0 .
%
%   Maximum drawdown is non-negative since it is the change from a peak to a
%   trough.
%
% References:
%	Malik Magdon-Ismail, Amir F. Atiya, Amrit Pratap, and Yaser S. Abu-Mostafa,
%	"On the Maximum Drawdown of a Brownian Motion", Journal of Applied
%	Probability, Volume 41, Number 1, March 2004, pp. 147-161. 
%
%	See also MAXDRAWDOWN

%	Copyright 2005 The MathWorks, Inc.
%	$Revision: 1.1.6.2 $ $Date: 2005/12/12 23:16:10 $

if nargin < 3 || isempty(Mu) || isempty(Sigma) || isempty(T)
	error('Finance:emaxdrawdown:MissingInputArg', ...
		'Missing required input arguments Mu, Sigma, or T.');
end

[n, i] = size(Mu);
if (n > 1) || (i > 1)
	error('Finance:emaxdrawdown:InvalidInputArg', ...
		'Argument Mu must be a scalar.');
end

[n, i] = size(Sigma);
if (n > 1) || (i > 1)
	error('Finance:emaxdrawdown:InvalidInputArg', ...
		'Argument Sigma must be a scalar.');
end

T = T(:);
[n, i] = size(T);
if i > 1
	error('Finance:emaxdrawdown:InvalidInputArg', ...
		'Argument T must be a vector.');
end
	
if ~isfinite(Mu) || ~isfinite(Sigma) || (Sigma < 0)
	error('Finance:emaxdrawdown:InvalidInputArg', ...
		'Invalid argument values. Requires finite Mu, Sigma >= 0, T >= 0.');
end

for i = 1:n
	if ~isfinite(T(i)) || (T(i) < 0)
		error('Finance:emaxdrawdown:InvalidInputArg', ...
			'Invalid argument values. Requires finite T >= 0.');
	end
end
 
if Sigma < eps
	if Mu >= 0.0
		EDD = zeros(n,1);
	else
		EDD = - Mu .* T;
	end
	return
end

EDD = zeros(n,1);

for i = 1:n
	if T(i) < eps
		EDD = zeros(n,1);
	else
		if abs(Mu) <= eps
			EDD(i) = (sqrt(pi/2) * Sigma) .* sqrt(T(i));
		elseif Mu > eps
			Alpha = Mu/(2.0 * Sigma * Sigma);
			EDD(i) = emaxddQp(Alpha * Mu * T(i))/Alpha;
		else
			Alpha = Mu/(2.0 * Sigma * Sigma);
			EDD(i) = - emaxddQn(Alpha * Mu * T(i))/Alpha;
		end
	end
end

function Q = emaxddQp(x)

A = [	0.0005;		0.001;		0.0015;		0.002;		0.0025;		0.005;
		0.0075;		0.01;		0.0125;		0.015;		0.0175;		0.02;
		0.0225;		0.025;		0.0275;		0.03;		0.0325;		0.035;
		0.0375;		0.04;		0.0425;		0.045;		0.0475;		0.05;
		0.055;		0.06;		0.065;		0.07;		0.075;		0.08;
		0.085;		0.09;		0.095;		0.1;		0.15;		0.2;
		0.25;		0.3;		0.35;		0.4;		0.45;		0.5;
		1.0;		1.5;		2.0;		2.5;		3.0;		3.5;
		4.0;		4.5;		5.0;		10.0;		15.0;		20.0;
		25.0;		30.0;		35.0;		40.0;		45.0;		50.0;
		100.0;		150.0;		200.0;		250.0;		300.0;		350.0;
		400.0;		450.0;		500.0;		1000.0;		1500.0;		2000.0;
		2500.0;		3000.0;		3500.0;		4000.0;		4500.0;		5000.0 ];

B = [	0.01969;	0.027694;	0.033789;	0.038896;	0.043372;	0.060721;
		0.073808;	0.084693;	0.094171;	0.102651;	0.110375;	0.117503;
		0.124142;	0.130374;	0.136259;	0.141842;	0.147162;	0.152249;
		0.157127;	0.161817;	0.166337;	0.170702;	0.174924;	0.179015;
		0.186842;	0.194248;	0.201287;	0.207999;	0.214421;	0.220581;
		0.226505;	0.232212;	0.237722;	0.24305;	0.288719;	0.325071;
		0.355581;	0.382016;	0.405415;	0.426452;	0.445588;	0.463159;
		0.588642;	0.668992;	0.72854;	0.775976;	0.815456;	0.849298;
		0.878933;	0.905305;	0.92907;	1.088998;	1.184918;	1.253794;
		1.307607;	1.351794;	1.389289;	1.42186;	1.450654;	1.476457;
		1.647113;	1.747485;	1.818873;	1.874323;	1.919671;	1.958037;
		1.991288;	2.02063;	2.046885;	2.219765;	2.320983;	2.392826;
		2.448562;	2.494109;	2.532622;	2.565985;	2.595416;	2.621743 ];

if x > 5000
	Q = 0.25*log(x) + 0.49088;
elseif x < 0.0005
	Q = 0.5*sqrt(pi*x);
else
	Q = interp1(A,B,x);
end

function Q = emaxddQn(x)

A = [	0.0005;		0.001;		0.0015;		0.002;		0.0025;		0.005;
		0.0075;		0.01;		0.0125;		0.015;		0.0175;		0.02;
		0.0225;		0.025;		0.0275;		0.03;		0.0325;		0.035;
		0.0375;		0.04;		0.0425;		0.045;		0.0475;		0.05;
		0.055;		0.06;		0.065;		0.07;		0.075;		0.08;
		0.085;		0.09;		0.095;		0.1;		0.15;		0.2;	
		0.25;		0.3;		0.35;		0.4;		0.45;		0.5;
		0.75;		1.0;		1.25;		1.5;		1.75;		2.0;
		2.25;		2.5;		2.75;		3.0;		3.25;		3.5;
		3.75;		4.0;		4.25;		4.5;		4.75;		5.0 ];

B = [	0.019965;	0.028394;	0.034874;	0.040369;	0.045256;	0.064633;
		0.079746;	0.092708;	0.104259;	0.114814;	0.124608;	0.133772;
		0.142429;	0.150739;	0.158565;	0.166229;	0.173756;	0.180793;
		0.187739;	0.194489;	0.201094;	0.207572;	0.213877;	0.220056;
		0.231797;	0.243374;	0.254585;	0.265472;	0.27607;	0.286406;
		0.296507;	0.306393;	0.316066;	0.325586;	0.413136;	0.491599;
		0.564333;	0.633007;	0.698849;	0.762455;	0.824484;	0.884593;
		1.17202;	1.44552;	1.70936;	1.97074;	2.22742;	2.48396;
		2.73676;	2.99094;	3.24354;	3.49252;	3.74294;	3.99519;
		4.24274;	4.49238;	4.73859;	4.99043;	5.24083;	5.49882 ];

if (x > 5.0)
	Q = x + 0.5;
elseif x < 0.0005
	Q = 0.5*sqrt(pi*x);
else
	Q = interp1(A,B,x);
end
